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ABSTRACT 

The  problem  of  radiowave  propagation  over  irregular  terrain  is  solved  by  using  the 

wide  angle  parabolic  equation.  The  terrain  is  characterized  by  its  height  profile  and  its 
ground  constants  (here  conductivity  a  — >  oo).  We  consider  horizontal  polarization  and 
treat  the  ground  as  perfectly  conducting  (PEC)  to  simplify  the  formulation.  This  thesis 

uses  a  piece-wise  conformal  transformation  to  flatten  the  irregular  terrain.  The  equations 

are  solved  by  the  split-step  Fourier  algorithm.  A  Harming  window  is  used  in  both  the 

spatial  and    wavenumber  domains  to  contain  the  computational  domain.  Effect  of  some 

numerical  parameters  such  as  the  horizontal  step  size,  and  height  of  the  computational 

domain  on  the  accuracy  of  the  solution  is  investigated.  The  numerical  results  are  compared 

with  available  results  for  some  typical  propagation  problems. 
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I.  INTRODUCTION 

A.  BACKGROUND 

Radiowave  propagation  over  irregular  terrain  and  in  the  presence  of  ducts  is  an 
extremely  important  topic  in  ground-to-ground  as  well  as  in  ground-to-air 
communications.  Similarly,  the  ability  to  predict  radiowave  propagation  over  irregular 
terrain  has  a  significant  impact  in  determining  target  detectability  in  a  radar  system.  The 
physics  of  propagation  is  affected  by  ever-changing  atmospheric  conditions  and  by 
complex  terrain  features  on  the  ground.  The  path  between  the  transmitter  and  receiver  is 
often  obstructed  by  natural  or  man-made  obstacles,  such  as  hills,  buildings,  atmospheric 
layers,  rain,  etc.  In  these  cases  waves  can  reach  the  receiver  along  more  than  one  path  and 
the  phenomenon  is  termed  as  multipath.  In  the  case  of  atmospheric  multipath  fading, 
super-refraction  or  sub-refraction  can  lead  to  abnormal  radiowave  propagation,  which  can 
result  in  high  gain  or  severe  loss  of  signal.  Reflection  multipath  fading,  which  is  due  to 
interference  between  the  direct  and  the  ground  reflected  waves,  depends  strongly  on  the 
terrain  geometry.  It  is  important  to  everyone  involved  with  communications  and  radar 
systems  to  predict  the  electromagnetic  fields  due  to  radiating  sources  in  the  troposphere 
taking  into  account  ducting  and  terrain  effects. 

Numerous  analytical  methods  are  available  for  predicting  electromagnetic  wave 
propagation,  such  as  geometric  optics,  physical  optics,  normal  mode  analysis  and 
combinations  of  the  above.  However,  a  complex  environment  complicates  the  application 
of  some  of  these  methods.  The  Parabolic  Equation  (PE)  method  has  emerged  as  an 
extremely  valuable  method  for  assessing  radiowave  propagation  in  the  lower  atmosphere 
in  the  presence  of  ducts  and  over  flat  terrain,  and  for  predicting  radiowave  propagation 
over  sloping  irregularities.  One  advantage  of  the  parabolic  Partial  Differential  Equation 
(PDE)  is  that  the  field  at  any  location  can  be  computed  in  terms  of  the  field  at  a  previous 
location.  However,  this  would  be  accurate  only  when  the  waves  propagate  predominantly 
in  the  forward  direction.  This  is  approximately  met  in  many  propagation  problems  and  is 
often  assumed. 


B.    OBJECTIVE 

In  this  thesis  we  adopt  a  wide  angle  parabolic  equation  due  to  Thomson-Chapman 
[Ref.  1  ]  to  predict  radiowave  propagation  over  an  irregular  terrain  and  in  the  presence  of 
ducts.  We  consider  horizontal  polarization  and  treat  the  ground  as  perfectly  conducting 
(PEC)  to  simplify  the  formulation.  In  practical  situations  this  assumption  will  have  minimal 
effect  for  frequencies  in  the  VHF  band  and  above.  The  key  feature  of  this  thesis  is  to  use 
piece-wise  conformal  transformations  to  map  the  irregular  terrain  into  a  flat  one  and  use 
known  techniques  available  for  flat  terrain.  It  is  in  the  transformed  space  rather  than  in  the 
physical  space  that  the  parabolic  approximation  made.  Chapter  II  presents  the  derivation 
of  the  governing  parabolic  equation  and  the  formulation  of  the  conformal  mapping. 
Chapter  III  details  the  generation  of  mesh  in  the  transformed  space  and  the  numerical 
procedure  for  solving  the  PDE.  The  accuracy  of  the  numerical  solution  is  examined  in 
Chapter  IV.  This  includes  a  study  on  the  effects  of  using  different  numerical  values  for 
various  important  parameters  (e.g.,  step  size,  maximum  height)  on  the  accuracy  of  the 
solution  and  validation  of  the  numerical  results  with  known  solutions  for  typical 
propagation  problems.  Recommendations  and  conclusions  are  presented  in  Chapter  V. 


II.  FORMULATION 

In  this  chapter  we  present  theory  on  the  governing  parabolic  equation  and 
conformal  mapping.  Figure  1  shows  a  horizontally  polarized  dipole  placed  over  an 
irregular,  PEC  terrain.  The  terrain  is  characterized  by  its  height  profile  and  its  ground 
constants  (here  the  conductivity  a  ->  «  ).  We  wish  to  solve  the  fields  at  a  point  over  the 
ground  in  the  presence  of  irregularities  and    ducts.  We  consider  a  2-dimensional  case 

where  the  source,  geometry,  refractivity  profile,  and  all  fields  are  y-invariant. 
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Figure  1 .  A  source  producing  fields  over  an  irregular  terrain. 


The  terrain  geometry  is  specified  by  means  of  discrete  points,  and  straight  line 
interpolation  is  assumed  between  the  given  data.  The  parabolic  equation  we  present  has  to 
be  finally  solved  numerically  using  either  finite  difference  or  Fourier  transform  approach. 
Both  of  these  methods  work  best  in  a  rectangular  domain  and  a  Cartesian  mesh.  Because 
the  lower  boundary  is  irregular,  we  cannot  adopt  a  Cartesian  mesh  without  discarding  or 
adding  points  in  going  from  one  range  step  to  the  next.  To  avoid  this  situation,  we  will  use 


a  coordinate  transformation  that  will  flatten  the  bottom  boundary.  Although  there  are 
multitudes  of  transformations  available  [Ref.  2],  we  have  decided  to  use  conformal 
mapping  in  our  case.  The  advantage  of  using  a  conformal  mapping  is  that  it  preserves  the 
Laplacian  operator  without  introducing  undesirable  cross-terms.  This  will  be  particularly 
advantageous  when  split-step  Fourier  transform  approach  is  used  to  solve  the  PE.  In  this 
thesis,  we  use  the  split-step  Fourier  transform  approach  to  solve  the  parabolic  equation 
due  to  its  superior  computational  perfomance.  Our  development  parallels  that  of  [Ref.  3]. 

A.  WIDE  ANGLE  PARABOLIC  PARTIAL  DIFFERENTIAL  EQUATION 

Figure  2  shows  the  coordinate  system  used  in  the  thesis. 
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earth's  surface 
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Figure  2.  Earth-centered  spherical  geometry. 

The  starting  point  for  our  formulation  is  the  earth  flattened  Helmholtz  equation  for 
the  scaled  field  V  [Ref  4]: 


V2xzV+k2o(\+2m)V=0 


(1) 


where   V(x,z)=  Jx  sin  0  E$(x,  z)   for   horizontal   polarization,   E^    being   the   azimuthal 
component  of  the  electric  field,  and  (r,6,(J))  are  the  usual  spherical  coordinates.  The  center 

of  the  coordinate  system  is  at  the  center  of  the  earth  and  angle  9  is  measured  from  the 


vertical  line  joining  the  center  and  transmitter  location.  The  transmitter  is  assumed  to  be 
located  at  <J)=0  (Figure  2).  Furthermore,  k0  =  co.ye0u,0  ,  is  the  free  space  wavenumber, 
x=rsin0  is  the  range  variable,  z  is  the  height  variable,  and  m(x,z)=(n-l+z/ae)  is  the  modified 
refractive  index  of  the  atmosphere  having  an  actual  refractive  index  n,  where  ae  is  the 
radius  of  the  earth.  In  this  chapter  we  include  only  the  final  details.  Complete  details  can 
be  found  in  [Ref  5]. 
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Figure  3.  Transformation  from  W-plane  to  ^-plane. 


Under  conformal  mapping,  [  Fig.  3],  equation  (1)  is  transformed  to: 


V^V+kllW  (Ql2(l+2m)F=0 


(2) 


On  expanding  the  Laplacian  operator,  equation  (2)  can  be  rewritten  as 


d2  V    d2  V  i 

■+-TT+k2o I W  (Ql2(l  +  2m) V=0  , 


2        Av,2 


dtf-      &r\ 


(3) 


where  I W  (Ql  is  the  Jacobian  of  the  transformation  and  is  given  by 


\W(Q\  = 


•    l-sin(^) 

'Li         .      .TCCJ 


(4) 


Substituting  V=Uelk°^ ,  equation  (3)  can  be  simplified  to: 


-^=ik0(Q-\)U 


(5) 


i  a2 

with  Q-l=i2      2 


1 


*"o+l»i$> 


+  m  ,  and  we  have  made  the  reasonable  assumption 


that    \W  (Q)\&  1.  For  the  transformation  considered  in  this  thesis,  this  is  true  everywhere 
except  near  points  P  and  Q  of  Figure  5. 

B.  CONFORMAL  MAPPING 

1.  Conformal  Mapping  Approach 

A    conformal    mapping   can   transform   an   irregular   polygonal    shape   into    a 
rectangular  one.  We  use  a  series  of  piece-wise  conformal  maps  to  transform  the  piecewise 
linear  terrain  to  a  flat  one  as  shown  in  Figure  4. 
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Figure  4.  Terrain  flattening. 
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For  example,  we  conformally  map  the  polygon  AABBA  shown  on  the  left  to  a 
corresponding  rectangle  shown  on  the  right.  Because  the  slope  is  different  from  one 
segment  to  the  other,  the  mapping  is  also  different.  Detailed  derivations  of  the  mapping 
are  given  in  [Ref.  5].  Here  we  present  only  the  final  results. 

2.  Analytical  Formulation  of  the  Mapping 

Consider  the  transformation  W=f(Q  between  the  W=r+iz  and  C=£,+ir|  complex 
planes  as  shown  in  Figure  5. 
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Figure  5.  Piece- wise  conformal  mapping. 

For  the  sake  of  simplicity  we  assume  that  the  points  R  and  S  are  at  infinity.  The 
mapping  is  given  by  the  Schwarz-Christoffel  formula  [Ref.  3]: 


w=Wp+%e™  J— i 


dX 


0  XV+2(l-}.)2 


i-v 


(6) 


1  n£> 

where  C2=2(sin(A^)+1)- 


With  reference  to  Figure  5,  we  define  the  slope  prameter  v  as 


*(a7)>  "2<V<2 


v=¥fc)>  -i<v<T-  (7) 


It  is  asumed  that  positive  angles  are  measured  in  the  counter  clockwise  direction.  Points 
on  the  line  PR  are  mapped  to  the  points  on  the  line  PR  through: 


z/=zy_1+^J(u;v)  (8) 


where  as  points  on  the  line  Q'S'  are  mapped  to  points  on  the  line  QS  through: 


Zr=ZJ+1fJ(u-v)  (9) 


2    wn      A    U 
where  u=tanh  (tt~)  =1 — ,       and 

v2Ar/      1-u 


u 


J(u;v)=j"  - 


,  x  x-v 


0   \2      (1-A,)2 

Accurate  formulas  for  computing  J(u;v)  are  given  in  [Ref.  5]. 


III.  SOLUTION  PROCEDURE 

In  the  previous  chapter  we  presented  the  underlying  equations  for  the  computation 
of  fields  over  irregular/flat  terrain  and  ducting  conditions.  In  section  A  of  the  present 
chapter,  we  start  with  the  generation  of  mesh  points  in  the  physical  space.  In  section  B,  we 
continue  with  the  numerical  implementation  using  the  split-step  Fourier  (SSF)  algorithm 
for  the  final  computation  of  the  fields. 

A.  GENERATION  OF  THE  TRANSFORMED  SPACE 

It  is  clear  from  equations  (8)  and  (9)  that  a  uniform  mesh  r|=kAr|  in  the 
transformed  domain  maps  to  different  distribution  of  points  on  the  vertical  lines  PR  and 

QS.  Furthermore,  the  distribution  of  points  is  dependent  on  the  slope  parameter  v. 

Therefore,  in  going  from  one  range  step  to  another,  the  distribution  of  points  on  the  right 

side  of  the  previous  step  will  not  coincide  with  the  distribution  of  points  on  the  left  side  of 

the  present  step.  This  is  illustrated  in  Figure  6  for  v  =  0.25,  Ar|  =  2  and  Ar  =  25  m,  where 


z=z  max=50  m 


z=0 


\z—25  m^r—25m^ 


Figure  6.  Distribution  of  points  about  a  slope  discontinuity. 


points  just  to  the  left  of  triangle  base  are  indicated  by  circles  and  points  to  the  right  of  it  by 
inverted  triangles.  Clearly  interpolation  or  extrapolation  of  data  is  needed  before  the  field 
values  calculated  at  one  range  can  be  used  as  input  for  the  next  step.  In  this  thesis,  we  use 
rational  approximation,  [Ref.  5]  to  do  the  desired  interpolation. 

B.  NUMERICAL  IMPLEMENTATION  USING  THE  SPLIT-STEP  FOURIER 
(SSF)  ALGORITHM 

Consider  the  parabolic  partial  differential  equation  given  in  equation  (5).  We 
would  like  to  implement  the  equations  by  using  the  SSF  algorithm.  Figure  7  shows  field  U 
on  a  given  vertical  line.  We  would  like  to  predict  the  field  U+  on  a  subsequent  line  spaced 
Ar  from  the  original  line. 
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77777777777777777? >  V=0 
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Figure  7.  Marching  in  range. 

To  solve  the  parabolic  equation,  the  following  Fourier  transform  pair  is  used  in  the 
transformed  domain: 


U &p)  =  F(U)  =  I  U&rti  e*™  di\ 


(10a) 
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D&rfl  =F\U)  =  £  Jb&p)  e-'^dp  (10b) 

The  boundary  condition  U(£,,0)=0  for  horizontal  polarization  implies  that 
U(£ ,-ti)=-U(£,ti),  T|>0,  and  U(^,-p)  =  -U(^,p),  p>0.  Using  equation  (10a)  in  equation 
(5)  and  making  use  of  standard  'splitting  of  operators'  [Ref.  2]  we  have: 

U+  =e~iAU-  (11) 

where        A=  .  .   The  field  in  the  spatial  domain  is  obtained  by  inverse 

h+Jfi-p2 

transformation  and  correcting  for  refractive  effects: 

U+  =  eik°mArF-l[e-iAF[U-]]  (12) 

In  computing  the  ordinary  Fourier  transform  pair,  we  use  a  N-point  FFT.  Let  us 
assume  that  the  various  quantities  are  band  limited  over  pmax  <  p  <  pmax  ,  and  that  the 
transform  is  evaluated  at  p=0,   Ap,   2Ap,    ....(N-l)Ap  where  Ap=27c/(NAz).   Positive 

wavenumbers  occur  at  p=Ap,  2Ap, ,(N/2-l)Ap  while  negative  wavenumbers  occur  at 

(N/2+1  )Ap,  (N/2+2) Ap, ,  (N-l)Ap.  The  value  (N/2)Ap  corresponds  to  both  +pmax  and 

"Pmax  We  use  a  Hanning  window  to  contain  the  computational  domain  at  the  top.  The 

Harming  window  is  given  by: 


1  ,0<«<^ 

sin   (— )     'T^"^7 


We  use  it  both  in  the  spatial  and  wavenumber  domains.  A  mirror  image  of  h(n)  about  n=0 
is  used  for  negative  wavenumbers.  Note  that  the  Hanning  window  forces  a  gradual  rolloff 
to  zero  over  the  last  quarter  of  the  domain.  With  the  use  of  the  Hanning  window  the 
solution  of  the  parabolic  equation  can  be  written: 

U+  =  elk°mbrF-x  [h^e-^Flh^U-]]  (14) 
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As  we  mentioned  at  the  begining  of  this  chapter,  because  consecutive  steps  are 
taken  in  different  transformed  spaces,  the  calculated  field  at  the  right  side  of  one  range 
step  cannot  be  directly  used,  in  general,  as  input  to  the  next  step.  For  example,  at  the 
junction  between  flat  and  sloping  terrain,  points  are  uniformly  distributed  over  a  flat 
surface,  but  unevenly  distributed  over  the  sloping  terrain.  As  mentioned  previously,  we  use 
the  rational  approximation  to  calculate  the  field  f(z)  at  a  desired  point  z,  in  terms  of  the 
fields  f  ,f+  and  f0  at  positions  z  ,  z+  and  z0  respectively  as  shown  in  Figure  8. 


Figure  8.  Field  interpolation. 


The  field  is  interpolated  as: 


f(z)  =  (f0+a1(z-zo))/(\  +  bl(z-zo)). 


(15) 


where      a ,  =  (f+s-  -fs+)/(f+  -/_) ,    b  i  =  (s=  -  s+)l(f+  -/-)    and     s-  =  <f0  -f-)/(zQ  -  z) , 


s+  =  (f+-f0)/(z+-z0). 
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IV.  RESULTS 

Having  introduced  the  basic  theoretical  background,  we  now  present  the  numerical 
results  using  the  developed  algorithm  for  several  practical  cases  of  radiowave  propagation. 
We  start  with  the  prediction  of  propagation  factor  (excess  loss  over  the  free  space  case) 
over  flat  earth,  with  and  without  a  duct.  We  compare  our  results  with  the  exact  solutions 
or  with  results  available  in  the  literature.  The  next  case  we  consider  is  that  of  propagation 
over  spherical  earth  in  a  non-refractive  atmosphere.  Finally,  we  present  our  solution  for 
propagation  over  a  knife-edge  shaped  obstacle.  We  use  this  terrain  irregularity  to 
comment  on  the  overall  accuracy  of  the  method  and  its  dependence  on  parameters,  such  as 
the  step  size  for  the  marching  and  the  maximum  height  chosen  to  truncate  the 
computational  domain. 

A.  PROPAGATION  OVER  A  PEC  PLANE 

Consider  a  transmitter  and  a  receiver  at  a  height  of  5  meters  above  flat  ground  and 
separated  by  1  km  ,  as  shown  in  Figure  9. 


J 

r~  i 

A 

r                                                                              ^ -\ 

B 

< 

=5  m 

J 

t_ 

< - 

|„  - ,  ,,,                                                                                   ;i 

k 


x=l  km 


H=5m 


Figure  9.  Propagation  between  two  stations  A  and  B,  both  of  which  are  stationed  over 
perfectly  conducting  ground. 
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In  our  numerical  approach,  the  propagation  factor,  PF,  is  obtained  from  the 
normalized  field,  V,  as  [Ref.  6] 

PF=\0\og(\V\2xXo)  (16) 

where  x  is  the  horizontal  from  the  transmitter.  A  positive  (negative)  value  of  propagation 
factor  implies  gain  (loss)  with  respect  to  propagation  in  free  space.  The  exact  solution  can 
be  obtained  from  image  theory  [Ref.  7,  p.  144].  The  expression  valid  in  the  far-field  and 
for  x  »  h,H  is  PF=\E/Ea\  =  2  sm(k0hcosQ),  where  cos  Q=H/jH2+x2  ,  h  is  the 
transmitter  height,  and  H  is  the  receiver  height.  The  transmitter  placed  at  x=0,  z=5,  where 

x,  z  are  in  meters.  The  frequency  of  operation  is  1  GHz.  We  choose  N=1024  ,  Az=Ar|= 
=0.15  m  and  Ax=Ar=5  m.  The  height  of  the  upper  boundary  is  zmax  =76.5  m.  The  electric 
field  is  determined  from  the  source  line  up  to  a  distance  of  1  km. 

The  variation  of  the  propagation  factor  versus  the  horizontal  distance  is  plotted  in 
Figure  10.  It  is  seen  that  there  is  a  very  good  agreement  with  the  image  theory  solution. 
For  the  first  50  m,  the  numerical  solution  differs  slightly  from  the  known  solution.  This  is 
due  to  the  breakdown  of  the  condition  x  »  h,  H  in  the  image  theory,  far-field  formula. 

B.  PROPAGATION  OVER  A  PEC  PLANE  IN  THE  PRESENCE  OF  DUCT 

Next  we  present  the  numerical  results  for  radiowave  propagation  over  a  PEC 
plane,  in  the  presence  of  a  tri-linear  duct  [Figure  1 1]  specified  by 


M(z) 


340  +  0.118z  0<z<135 

499.03  -1.06z      135<z<150  (17) 

322.33  +  0.118z  z>150 


where  z  is  in  meters  and  M  =  m  x  106.  This  is  a  surface  elevated  duct  with  AM=16  and 
Ah=150  m,  where  AM  and  Ah  are  defined  at  Figure  11.  The  minimum  frequency  (cut  off 
frequency)  required  for  propagation  in  the  presence  of  this  duct  is  217  MHz.  In  our  case 
the  frequency  of  operation  of  3  GHz  is  well  above  this  cut  off  frequency.  The  height  of  the 
transmitter  is  30  m  and  the  maximum  height  is  chosen  to  be  z       =512  m.   We  choose 

*-*  max 
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Figure  10.  Propagation  factor  versus  the  horizontal  distance  for  a  source,  placed  at  a  height  of 
5  m.  Flat  earth,  f =1  GHz  ,  Ax=5  m  and  H=h=5  m. 
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Figure  1  l.Tri-linear  surface  elevated  duct. 
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Figure  12.  Propagation  factor  versus  the  receiver  height  at  a  range  of  40  km  in  the  presence  of  a  tri-linear 
duct.  Flat  earth,  transmitter  height  =30  m,  f  =3  GHz. 
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N=512,  Dz=2  m  and  Dx=200  m.  Figure  12  shows  propagation  factor  versus  the  receiver 
height  at  a  distance  of  40  km.  A  favorable  agreement  with  [Ref.6,  Figure  3]  is  obtained. 

C.  PROPAGATION  OVER  SPHERICAL  EARTH 

Figure  13  shows  a  transmitter  and  a  receiver  located  beyond  the  horizon  over 
spherical  earth.  We  consider  transmitter  and  receiver  antennas  at  heights  of  10  m  each. 
Distance  to  the  horizon,  dT  ,  from  each  station  is  JlhRe  ,  where  Re  is  the  effective  radius 
of  the  earth  equal  to  6375  km.  In  our  case  a^  =11.3  km.  The  situation  described  above  is 
equivalent  to  the  one  where  radiowaves  propagate  over  flat  earth  but  in  the  presence  a 
modified  atmosphere  described  by  modified  index 


m(z)  =  (»-!+  zlae) 


(18) 


where  ae  is  the  radius  of  the  earth,  and  we  choose  the  refractive  index,  n,  equal  to  1 .  We 
use  the  equivalent  model  in  our  numerical  calculations. 


n=l 


m=z/a. 


transmitter 


Original  Problem 


receiver 


transmitter 


h=10m 


receiver 
y 

H=10m 


::::'::•:::•:•'••••••• """  '"■'". :'.':. :::;:":::.'. :.'.-*! 


d=  10  km 


Equivalent  Problem 


Figure  13.  Propagation  over  spherical  earth. 

We  consider  a  frequency  of  300  MFfz  and  treat  the  earth  as  a  PEC.  The  exact  solution  to 
the  problem  can  be  obtained  by  using  the  mode  theory  [Ref.  8].  Our  case  pertains  to  the 
case  of  low  antennas  where  the  heights  are  less  than  the  critical  height  hc ,  defined  by 
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hc  -  30Xo    ,  A.    in  m  . 


(19) 


In  our  case  hc=30  m  at  f=300  MHz.  The  distance  d,  must  satisfy:  X,0d  >  2  h  H,  meaning 
that  the  two  stations  must  be  at  least  200  m  apart.  We  use  N=1204,  Ax=50  m  and  Az=0.5 
m.  This,  places  the  maximum  height  at  zmax=255.5  m.  We  solve  for  the  propagation  factor 
as  a  function  of  the  horizontal  distance  using  our  model  and  compare  the  result  with  the 
exact  solution  [Ref.  8] 


PF   =  \EIE0\=lJ7^\  S 


5+2t„ 


fnW„(H)\ 


(20) 


where  5  =  [k0Re]^(zrc-  1),  src=er-iar   where  er  is  the  relative  permittivity  and  ar  is  the 

i 
relative  conductivity,  C,  =  sd,  s  =  [k0IR%\  3  ,xn  =  complex  numbers  which  depend  on  8,  and 

fn  (h),  fn(H)  are  height  gain  functions  for  the  n-th  mode  (term).  For  low  antennas  fn  (x)  is 
independent  of  n  and  is  given  by   fix)  =  1  +ik0xjerc-  1  -  The  contants  tn  are  given  as  in 

Table  1: 


n 

\ 

1 

1.856e"i7t/3 

2 

3.245e"i7t/3 

3 

4.382ei,l/3 

>3 

0.5[37r(n+0.25)f3e"i7t/3 

Table  1.  Mode  numbers  xn  for  8=00  ,(PEC). 

In  Figure  14  we  plot  the  propagation  factor  versus  the  horizontal  distance  d  for  a 
range  up  to  10  km.  There  is  a  very  good  agreement  with  the  modal  solution,  beyond  the 
first  400  m.  It  was  already  mentioned  that  the  distance  between  the  two  stations  has  to  be 
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Figure  14.  Propagation  factor  versus  horizontal  distance  over  spherical  earth,  f  =300  MHz,  Ax=50  m. 
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greater  than  200  m  in  the  modal  solution,  which  explains  the  disagreement  for  the  first  few 
hundred  meters. 

D.  PROPAGATION  OVER  PEC  KNIFE-EDGE 

The  final  example  we  consider  is  that  of  propagation  over  a  knife  edge  on  a 
conducting  plane.  Consider  a  perfectly  conducting  knife-edge  of  height  25  m,  as  shown  in 
Figure  15. 


h=3  m 


r — "• 


1 "■,,,,,,,,, ,,,,,,,„,,.,,,,,,,,,,,,,,,,,,,,,,,,,,„,,„,     >m1      ,    ,,.,,,,, ,,,.,,..,..y. mi    ,„,, ,,,,, , 


25  m 


receiver 
height 


B 


500  m 


<25  mJ<2g  *£ 


500  m 


Figure  15.  Perfectly  conducting  knife  edge  between  the  transmitter  at  A  and  the  receiver 
B,  both  of  which  are  on  perfectly  conducting  ground. 


We  compare  our  numerical  solution  with  a  four  ray  model  of  knife  edge  diffraction 
theory  [Ref.  9].  The  two  stations  are  at  distances  dx  and  d2  from  the  plane  of  the  edge.  We 
choose  d,=d2  =500  m.  In  the  parabolic  equation  computations,  we  replace  the  knife  edge 
with  a  triangular  hill  having  a  base  of  50  m  (Figure  15).  The  source  is  placed  at  a  height 
of  3  m  and  we  would  like  to  determine  the  propagation  factor  as  a  function  of  receiver 
height.  We  choose  a  maximum  height  of  255.5  m,  and  the  other  parameters  as:  N=1024, 
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Az=0.5  and  Ax=l  m.  The  frequency  of  operation  is  300  MHz.  The  numerical  results  are 
plotted  in  Figure  16,  along  with  the  four-ray  solution,  [Ref.  9].  There  is  an  excellent 
agreement  between  the  two  plots. 

E.  PARAMETRIC  STUDY  OF  THE  SSF  ALGORITHM  USING  CONFORMAL 
MAPPING 

We  use  the  knife  edge  irregularity  to  show  the  dependence  of  the  solution  on  the 
following  parameters:  the  horizontal  step  size  and  the  maximum  height  chosen.  We  choose 
three  different  heights,  each  of  which  corresponds  to  the  following  values  of  N  (with  fixed 
vertical  step  size  Az=0.5  m):  N=256,  512  and  1024.  The  values  of  height  are:  zmax=63.5, 
127.5  and  255.5  meters  respectively.  For  each  of  these  heights,  we  create  three  sets  of 
plots,  for  three  different  values  of  the  horizontal  step  sizes:  Ax=  5,  2.5  and  0.75  meters.In 
Figures  17,  20  and  23  (Ax=5  m),  we  discuss  the  effect  of  the  maximum  height  chosen  for 
fixed  Ax  and  it  is  seen  that  zmax  =76.5  m  (Figure  17)  is  certainly  not  high  enough.  Only  if 
we  use  higher  a  computational  domain  (Figures  20  and  23)  do  we  obtain  a  solution  which 
has  resemblance  with  the  four-ray  solution.  In  the  second  set  of  Figures  for  a  fixed  Ax=2.5 
m,  (Figures  18,  21  and  24),  again  we  conclude  that  the  higher  the  zmax  the  better  the 
agreement  between  the  numerical  and  the  four  ray  solution.  In  Figures  21  and  24  (zmax 
equal  to  125.5  and  255.5,  respectively)  we  start  to  obtain  a  solution  similar  to  the  four-ray 
theory  solution,  but  not  yet  an  accurate  one.  Finally  in  Figures  19,  22  and  25,  for  the  fixed 
horizontal  step  size  Ax=0.75  m,  we  see  that  the  solution  with  zmax  of  125.5  m  or  255.5  m 
(Figures  22  and  25)  produces  very  good  agreement. 

Next  we  discuss  the  effect  of  the  horizontal  step  size  Ax  for  a  fixed  vertical  step 
size  Az  and  zmax  (N).  In  Figures  17  through  19  the  maximum  height  is  76.5  m  and  Az=0.5 
m.  We  see  that  even  for  very  small  Ax,  an  acceptable  solution  is  not  obtained.  In  this  case 
the  solution  is  dominated  by  the  height  parameter  (not  convergent  yet).  In  Figures  20 
through  22,  we  use  zmax=125.5  m  and  Az=0.5  m.  It  is  seen  that  as  we  decrease  the 
horizontal  step  size  a  better  agreement  with  the  four-ray  solution  is  achieved  (Figure  22). 
Figures  23  through  25  present  the  results  for  a  fixed  Az=0.5  m  and  zmax  =255.5  m.  Again 
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Figure  16.  Propagation  factor  versus  the  receiver  height  at  a  range  of  1  km.  Smooth  earth,  the 
source  placed  at  a  distance  of  500  m  from  the  edge.  Az=0.5  m,  Ax=0.75  m  and  zmax  =255.5  m. 
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Figure  17.  Propagation  factor,  PF,  versus  the  receiver  height.  Az=  0.5  m,  Ax=5  m  and  zmax  =63.5  m. 
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Figure  18  Propagation  factor  versus  the  receiver  height.  Az=  0.5  m,  Ax=2.5  m  and  zmax  =63.5  m. 


25 


70 


60 


50 


*CD 

X 


40 


CD 

I  30 

o 

CD 


20 


10 


0 
-30 


Numerical  Solution 
Four  Ray  Theory 


-25 


-15  -10  -5 

Propagation  Factor  (dB) 


Figure  19.  Propagation  factor  versus  the  receiver  height.  Az=  0.5  m,  Ax=0.75  m  and  zmax  =63.5  m. 
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Figure  20.  Propagation  factor  versus  the  receiver  height  Az=  0.5  m,  Ax=5  m  and  zmax  =127.5  m. 
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Figure  21.  Propagation  factor  versus  the  receiver  height.  Az=  0.5  m,  Ax=2.5  m  and  z      =127.5  m. 
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Figure  22.  Propagation  factor  versus  the  receiver  height.  Az=  0.5  m,  Ax=0.75  m  and  zmax  =127.5  m. 
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Figure  23.  Propagation  factor  versus  the  receiver  height.  Az=  0.5  m,  Ax=5  m  and  zmax  =255.5  m. 
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Figure  24.  Propagation  factor  versus  the  receiver  height.  Az=  0.5  m,  Ax=2.5  m  and  zmax  =255.5  m. 
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Figure  25  Propagation  factor  versus  the  receiver  height.  Az=  0.5  m,  Ax=0.75  m  and  zmax  =255.5  m. 
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the  effect  of  Ax  is  easily  seen  and  Figure  25  shows  a  numerical  solution  with  an  excellent 
agreement  with  the  four-ray  theory  solution. 
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V.  CONCLUSIONS 

In  this  thesis,  a  numerically  efficient  method  to  model  radiowave  propagation  over 
irregular  terrain  and  in  a  stratified  atmosphere,  using  the  parabolic  equation,  was 
implemented  and  tested.  The  parabolic  equation  method  is  a  full-wave  method,  and 
aspects  such  as  forward  reflection,  refraction,  diffraction,  and  surface  wave  propagation 
are  included.  However,  it  ignores  back-scattering.  Since  it  ignores  back  propagation,  it 
allows  for  a  rapid  solution  of  the  fields  using  marching  techniques  along  the  propagation 
path  starting  from  the  source.  It  offers  an  advantage  compared  to  the  ray  method  in  that  it 
is  valid  in  the  shadow  region  where  the  latter  method  completely  breaks  down. 
Furthermore,  it  is  the  most  practical  method  for  predicting  propagation  over  long  ranges 
(thousands  of  wavelengths). 

After  presenting  the  wide  angle  parabolic  equation,  we  presented  a  piece-wise 
conformal  transformation  to  map  the  irregular  terrain  into  a  flat  one,  so  that  known 
techniques  available  for  flat  terrain  can  be  used.  The  split-step  Fourier  (SSF)  algorithm 
was  used  in  the  computational  domain  to  solve  the  parabolic  PDE.  A  Harming  window 
was  used  both  in  the  spatial  and  wavenumber  domains.  Although  very  small  step  sizes 
were  used  for  marching  along  ranges  of  several  kilometers,  the  method  is  very  time 
efficient.  Horizontal  polarization  was  considered  and  the  ground  was  treated  as  a  perfect 
electric  conductor  (PEC)  to  simplify  the  formulations.  This  assumption  will  have  minimal 
effect  in  practical  situations  and  for  frequencies  in  the  VHF  band  and  above. 

Numerical  results  to  predict  radiowave  propagation  for  several  practical  cases 
were  computed  and  validated  using  the  results  available  in  the  literature.  These  included 
both  propagation  over  flat  terrain  in  a  refractive  and  a  non-refractive  atmosphere, 
propagation  over  spherical  earth,  and  propagation  over  a  PEC  knife  edge  placed  on  a 
perfectly  conducting  plane.  Excellent  agreement  was  observed  and  demonstrated  for  all 
cases.  Also  included,  was  a  parametric  study  of  the  algorithm  used,  to  show  the  effect  of 
the  horizontal  step  size  Ax  and  the  maximum  height  zmax  ,for  a  fixed  vertical  step  size  Az. 
The  need  for  larger  heights  and  smaller  horizontal  steps  was  demonstrated. 

The  model  presented  in  this  thesis  can  be  applied  to  communication  problems  or 
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for  target  detection  where  prediction  of  propagation  of  HF/VHF  signals  over 
inhomogeneous  and  irregular  terrain  is  highly  desirable.  Horizontal  polarization  applies 
mainly  for  radar  systems.  For  this  case,  signals  received  from  a  target  depend  on  the  direct 
wave  and  the  ground  reflected  wave.  The  latter  ray's  path  depends  on  the  terrain 
roughness  and  ground  constants.  It  is  also  important  to  assess  the  atmospheric  conditions. 
The  numerical  model  developed  in  this  thesis  provides  a  useful  tool  for  predicting  the 
performance  of  a  radar  or  communication  system  before  its  final  development. 
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